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Abstract. The MUSIC algorithm, and its extension for imaging sparse extended objects, with 
noisy data is analyzed by compressed sensing (CS) techniques. A thresholding rule is developed 
to augment the standard MUSIC algorithm. The notion of restricted isometry property (RIP) 
and an upper bound on the restricted isometry constant (RIC) are employed to establish sufficient 
conditions for the exact localization by MUSIC with or without noise. 

In the noiseless case, the sufficient condition gives an upper bound on the numbers of random 
sampling and incident directions necessary for exact localization. In the noisy case, the sufficient 
condition assumes additionally an upper bound for the noise-to-object ratio in terms of the RIC 
and the dynamic range of objects. This bound points to the superresolution capability of the 
MUSIC algorithm. Rigorous comparison of performance between MUSIC and the CS minimization 
principle. Basis Pursuit Denoising (BPDN), is given. 

In general, the MUSIC algorithm guarantees to recover, with high probability, s scatterers with 
n = 0{s^) random sampling and incident directions and sufficiently high frequency. 

For the favorable imaging geometry where the scatterers are distributed on a transverse plane 
MUSIC guarantees to recover, with high probability, s scatterers with a median frequency and 
n = 0{s) random sampling/incident directions. 

Moreover, for the problems of spectral estimation and source localizations both BPDN and 
MUSIC guarantee, with high probability, to identify exactly the frequencies of random signals 
with the number n = 0{s) of sampling times. However, in the absence of abundant realizations 
of signals, BPDN is the preferred method for spectral estimation. Indeed, BPDN can identify 
the frequencies approximately with just one realization of signals with the recovery error at worst 
linearly proportional to the noise level. 

Numerical results confirm that BPDN outperforms MUSIC in the well-resolved case while the 
opposite is true for the under-resolved case, giving abundant evidences for the superresolution 
capability of the MUSIC algorithm. 

Another advantage of MUSIC over BPDN is the former's flexibility with grid spacing and guar- 
antee of approximate localization of sufficiently separated objects in an arbitrarily refined grid. 
The localization error is bounded from above by 0(\a) for general configurations and by 0{\) for 
objects distributed in a transverse plane. 



1. Introduction 

The MUSIC (standing for MUltiple-Signal-Classification) algorithm is a well-known method in 
signal processing for estimating the individual frequencies of multiple time-harmonic signals [Sj [21] . 
Mathematically, MUSIC is essentially a method of characterizing the range of the covariance matrix 
of the signals (see Section [6] for details). 

MUSIC was originally developed to estimate the direction of arrival for source localization [19j . 
Later, the MUSIC algorithm is extended to imaging of point scatterers [BJ. A proof of a sufficient 
condition for the exact recovery of the object support in the noiseless case is given in [16] (see 
also [15j) which is reproduced in Proposition [l] below. The performance guarantee is general but 
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qualitative in nature. Neither does it take into account the effect of noise which is important for 
assessing the superresolution effect. 

The main purpose of this paper is to give a quantitative performance evaluation for the MUSIC 
algorithm in terms of how many data are needed and how they may be collected in order to exactly 
recover the locations of given (large) number of objects, be they sources, scatterers or frequencies 
as well as how much noise the MUSIC algorithm can tolerate. Our approach is based on recent 
advances in compressed sensing theory ( [H [21 [18] and references therein) and applications to imaging 
( [HI [ini [Hj and references therein) . 

A main result for localizing scatterers obtained in the present paper has the following flavor 
(more details later): Let ^max and ^min be, respectively, the strengths of the strongest and weakest 
(nonzero) scatterers, Sf- the (generalized) restricted isometry constants (RIC) of order s and e the 
level of noise in the data. If the noise-to-scatterer ratio (NSR) obeys the upper bound 



(1) ^ < ^/(l + <5+)2^ + (l-5.-)2A-(l + 5 



where 

A 1 1 ' 



2 2 7V2rs + 1 



and (defined in (19)) is a measure of the independence of the column vectors outside the 
object support from the range of the data matrix then the MUSIC imaging function with the 
thresholding rule 

recovers exactly the locations of the s scatterers (cf. Theorem [2| Section 3). Compressed sensing 
theory comes into play in addressing the dependence of RIC on the frequency, the number and 
distribution of random sampling directions (or sensors), the number of scatterers and the inter- 
scatterer distances. 

In the super-resolution regime, the 5j tends to 1 and Ts tends to zero, rendering the right hand 
side of ([l]) approximately 

(3) '^-/-''^ 

2(1 + 5J)emax/e mm 

where ^max/'Cmin IS the dynamic range of the scatterers. For a NSR smaller than ^ the s scatterers 
can still be perfectly localized by the MUSIC algorithm with the thresholding rule ^ where the 
threshold is approximately 

» 1. 



(l-C+i)' V2 + 5s^ 

Previous observation ^17j and our numerical results (Section [s]) lend support to this superresolution 
effect of the MUSIC algorithm. 

First let us review the inverse scattering problem and the MUSIC imaging method. 
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Figure 1. Far- field imaging geometry 



1.1. Inverse scattering. Consider the scattering of the incident plane wave 

(4) n'(r) = e*^""-^ 

by the variable refractive index n^(r) = 1 + ^(r) where d is the incident direction. The scattered 
field satisfies the Lippmann-Schwinger equation 

(5) u%r) = 00^1 ^{r'){u\r')+u'{r'))G{r,r')dr', d=2,3 

where G(r,r') is the Green function of the operator —(A + w^) [16j. We assume that the wave 
speed is unity and hence the frequency equals the wavenumber a;. 
The scattered field has the far-field asymptotic 

(6) n^(r) = ^-^^(A(?,d) + 0(|rri)), f = r/|r|, 
where the scattering amplitude A is determined by the formula 

(7) A(r,d) = ^ / dr'ar'Mry-^'''^'-'. 

In the Born regime, the total field u on the right hand side of ([T]) can be replaced by the incident 
field u\ 

The main objective of inverse scattering then is to reconstruct the medium inhomogeneities ^ 
from the knowledge about the scattering amplitude A{r,d). 

Next we recall the MUSIC algorithm as applied to localization of point scatterers. 

1.2. MUSIC for point scatterers. Let S = {rj : j = 1, ...,s} be the locations of the scatterers. 
Let 7^ 0,j = l,...,,s be the strength of the scatterers. We will make the Born approximation 
first and discuss how to lift this restriction at the end of the section (Remark [2]) . For the discrete 
medium the scattering amplitude becomes the finite sum 

(8) ^(r,d) = ^EGV(r 
under the Born approximation. 



Let d/,/ = l,...,m and Sfc,A; = l,...,n be, respectively, the incident and sampling directions. 
For each incident field di,l = l,...,m, the scattering amplitude is measured in all n directions 
Sk,k = 1, ...,n. The whole measurement data consist of the scattering amplitudes for all pairs of 
(di,Sfc). 

Define the data matrix Y = (Yfc,/) G C"^™ as 



(9) 



Ykj A{sk,di), k=l,...,n, l = l,...,m 



where we keep open the option of normalizing Y in order to simplify the set-up. The data matrix 
is related to the object matrix 

X = diag(e,) G j = l,...s 

by the measurement matrices $ and ^ as 
(10) Y = *X** 

where $ and ^ are, respectively, 



(11) 
(12) 



1 



n 



Id 



n 



after proper normalization. Both (11) and (12) are normalized to have columns of unit 2-norm. 



We extend the formulation (10)-(12) to the case of sparse extended objects in Appendix [A} 

Note that both $ and ^ are unknown and (10) can be inverted only after the locations of 



scatterers are determined. This is what the MUSIC algorithm is designed to accomplish. 

The standard version of MUSIC algorithm deals with the case of n = m and Sfc = dfc. A; 
as stated in the following result. 



l,...,n 



Proposition 1. [151 [T6] Let {s^ = dk,k G N} 6e a countable set of directions such that any 
analytic function on the unit sphere that vanishes in Sfc,VA; G N vanishes identically. Let /C C 
he a compact subset containing S. Then there exists no such that for any n > no the following 
characterization holds for every r £ fC: 

1 



n 



(13) r G 5 if and only if (j)y. 
Moreover, the ranges of $ and Y coincide. 



-iujSn-r\T 



G Ran(*). 



Remark 1. As a consequence, r £ S if and only if'P(pr = where V is the orthogonal projection 
onto the null space of^* (Fredholm alternative). And the locations of the scatterers can he identified 
by the singularities of the imaging function 



1 



(14) J(r) 

Moreover, once the locations are exactly recovered, then both $ and ^ are known explicitly and 
the strength = l,...,s of scatterers can he determined by inverting the linear equation (10) 
which is an over- determined system. 
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Remark 2. The assumptions of Proposition\^can be relaxed: instead of = dfc,VA:, it suffices to 
have ^ S (j^mxs y;/^^^,/^ /^^g rank s. 

In light of this observation, it is also straightforward to extend the performance guarantee for the 
Born scattering case to the multiple-scattering case. In the latter case, ^ consists of entries which 
are the total field evaluated at rj for the incident direction di, i.e. 

(15) ^i^j = u*{vj-di). 

What is really needed is that ^ G £m^s rank s since then Z = X^'* and X share the same 
support (see more on this in Section^. Generically this is true for sufficiently large m as we will 
show below. 

Define the incident and full field vectors at the locations of the scatterers: 

U\di) = (n'(ri;dO,...,n'(r.;dOf 
U{di) = {u{rr,di),...,u{rs:di)f €a. 

Denote 

(16) G = [(l-<5,j)G(rj,r,)]Ge^^ 

The discrete version of the Lippmann-Schwinger equation (i.e. the Foldy-Lax equation) can be 
written as 

(17) U{di) = U\di)+u;^GXU{di), l = l,...,m. 



The 5ij terms in (16) represent the singular self-energy terms of point scatterers and should be 
removed for self consistency. 

Denote U' = [U'{di), ...,U'{dm)] G C^""" and U = [[/(di ),..., C/(d„)] G C^^™. Suppose that 



uj is not an eigenvalue o/GX. Then we can invert eq. (11) to obtain 



U = (I- w2gx)-^u\ 

Hence U has rank s if does. Indeed, for sufficiently high frequency uj and m randomly se- 
lected incident directions with sufficiently large ratio ^Jrnjs, U' has rank s with high probability 
(Propositions^ and^ below). 

For some special imaging geometry it is possible to reduce the number of incident and sampling 
directions to 0{s) (Section^. 

1.3. Outline. Proposition [T] says that if the number of sampUng directions is sufficiently large then 
the locations of the s scatterers can be identified by the s singularities of J. However, the condition 
is only qualitative in the sense that an estimate for the threshold no is not given. It would be of 
obvious interest to know, e.g. how no scales with s when s is large and when the conventional 
wisdom (no = s + 1) derived from counting dimensions is true. Also, how much noise can the 
MUSIC algorithm tolerate? 

However, unless additional constraints are imposed on the measurement scheme (the frequency, 
the incident and sampling directions etc), it is unlikely to make progress toward obtaining an 
useful estimate which is the objective of the present study. In [6] a geometric constraint on the 
configuration of sensors and objects has been pointed out for exact recovery in the absence of noise. 
Moreover, it seems possible that a non-vanishing portion of s randomly distributed scatterers may 
not be exactly recovered in the presence of machine error no matter how large n is (Figure [5| middle 
panel, and Figure [7| right panel). 
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Let us briefly sketch our approach and results: We shall discretize the problem by using a finite 
grid for the computation domain fC and put the problem in a probabilistic setting by using random 
sampling directions. Moreover, we consider noisy data and aim for a result for stable recovery by 
MUSIC. For the case of well-resolved grids, we show by using the compressed sensing techniques 
that for the NSR obeying and with high probability, uq = 0{s'^) for general configuration of s 
objects and uq = 0{s) for objects distributed on a transverse plane. For the case of under-resolved 
grids, we seek sufficient conditions for approximate, instead of exact, localization of objects and we 
show that for sufficiently small NSR and with high probability, the localization error is 0{Xs) with 
n = 0(5^) for a general object configuration and the localization error is 0{X) with n = 0{s) for 
objects distributed in a transverse plane. 

Our plan for the rest of the paper is to first give a sensitivity analysis for MUSIC and derive 
the condition for exact recovery with noisy data under which the MUSIC algorithm based on the 
perturbed data matrix can still recover exactly the object support (Section[2]). Next, we review the 
basic notion of compressed sensing (CS) theory and show how it naturally lends itself to a proof 
of exact localization by MUSIC (Section [s]). We show that with generic, random sampling and 
sufficiently high frequency the MUSIC algorithm can, with high probability, recover s scatterers 
with n = O(s^) sampling and incident directions (Corollary [2]). Then we consider a favorable 
imaging geometry where the scatterers are distributed on a transverse plane (Section |4]). We show 
that with median frequency the MUSIC algorithm can recover, with high probability, s scatterers 
with n = 0{s) sampling and incident directions (Corollary [3]). Next we analyze the performance 
guarantee of the compressed sensing principle. Basis Pursuit Denoising (BPDN) (Section [5]) and 
show that in the generic situation BPDN with sufficiently high frequency can recover s scatterers 
with n = C(s^) sampling directions and just one incident wave (Remark 10) while for the favorable 
geometry of planar objects BPDN can recover s scatterers with n = 0{s) sampling directions 
and one incident wave (Remark [9| . In Section [6] we return to the original applications of MUSIC 
and perform the compressed sensing analysis of the performance of MUSIC as applied to spectral 
estimation and source localization. We show that the MUSIC algorithm can, with high probability, 
identify exactly the frequencies of random signals with the number n = 0{s) of sampling times 
(Corollary [g] and Remark We discuss MUSIC in the setting with an arbitrarily fine grid and 
give error bounds in Section [7j Numerical tests are given in Section [8] where the superresolution 
capability of MUSIC and the noise sensitivity are studied. We conclude in Section |9j We give an 
extension of the MUSIC algorithm to the case of extended objects in Appendix |^ and a proof of 
performance guarantee for BPDN in Appendix |B| 



2. Sensitivity analysis 

For quantitative performance analysis of the MUSIC algorithm, we will work with the discrete 
setting and assume that /C is a discrete set of A^, typically large, number of points, i.e. the com- 
putation grid. The discrete setting appears naturally in applying MUSIC to imaging of extended 
scatterers (see Appendix [a|). Moreover, we consider the extension ^ of $ which includes not only 
the columns (t>rj,j = 1; representing the locations of the objects but also the columns repre- 
senting all the points in /C. Hence ^ G C"^^ and as usual $ is normalized so that the columns 
have unit 2-norm. The ordering of the columns of $ is not important for our purpose as long as 
they correspond to the points in /C in a well-defined manner. G (^nxN jg gij^tilarly de&ied. Also 
the extension X G C^^^ of X is defined by filling in zeros in all the entries outside the object 
support. 

6 



In terms of these notations, we can write 

N 

(18) Y = *X** = ^ ® 

i=i 

By a slight abuse of notation, we shall use S to denote the locations of objects in the physical 
domain as well as the corresponding index set. Likewise 5^ denotes the complement set of S in 
the computation grid fC as well as the total index set {1, A^}. In the same vein, $5 denotes the 
column submatrix of $ restricted to the index set S. Hence $5 = $ and = ^• 

First, let us reformulate the condition (13) for exact recovery as follows. 

Note that 

is the orthogonal projection of (pr onto the range of $ where is the pseudo- inverse of Hence 

V^r = (I - **^)</'r. 

If (pr for r G S'^ is independent of the columns of then 



and vice versa. Therefore ( 13 ) is equivalent to 



(19) = min ll^rlls ^ll^'/'rlb = ,/l -max||,/.r||-2,^^**t,^r > 0. 

reS'^ V rSiS'^ 

The number Ts gives a measure of how "independent" (j)z is from the range of $ uniformly in 
r G S^. 

Now we give a sensitivity analysis for MUSIC with respect to perturbation in the data matrix 
Y in terms of and other parameters. We want to show what else is needed, in addition to (13), 
to guarantee exact recovery of the support of scatterers when the data matrix is perturbed. 

The general data matrix considered in this paper has the form Y^ = Y + E where Y = $X'J'* G 
C"^"*, m > s, the number of objects. Set Z = X^'* G C**^"^ such that Y = $Z where Z is assumed 
to have rank s. 

We shall treat Z as the new object matrix and consider perturbed data matrices of the form 

(20) Y^ = *Z + E. 

Note that the locations of objects represented by Z are identical to those represented by X = 
diag(^j). 
Set 

ye ^ Y^Y^* = y + £ 

where 

(21) y = *zz*** G c"''" 

(22) £ = EZ*** + *ZE* + EE* G C"""" 

are both self-adjoint. Note that the range of y is the same as the range of $ and under the 
assumption of (13) equals to the span of {(pr : r G S}. 

Let {vj : j = 1, s} and {v^ : j = s + 1, n}, respectively, be the set of orthonormal bases for 
the range and null space of y. Let Qi G C"^^ and Q2 G respectively, be the matrices 

whose columns are exactly {vj : j = 1, s} and {vj : j = s + 1, n}. Let Q = [Qi, Q2] G C"^". 



Let o"! > (72 > • • • > o"„ be the singular values of y. Denote the smallest nonzero singular value 
of y by CTmin and set cJmax = ci- If Y has rank s, then Umin = Cs- We partition Q*<£'Q as follows: 



(23) q*£:q 



£11 £12 
£21 £22 



where £u e C^^^^l2 e c^^("-^),£:2i e c("-^)^^£:22 G c("-^)^("-^). 

The following is a slight recasting of a general result of matrix perturbation theory [20] . 
Proposition 2. (Theorem 2.7, Chap. V, ^) If 

(24) ^^^t^^f't 11 - < ^ 

f^mm — ||ill||2 — ||i22||2 ^ 

t/ien t/iere exisi F G c("-'5)><'' ^i/i 

(25) IIFII2 < ^ff' „^ „ 

Cmin — ||t^ll||2 — ||t^22||2 

such that the columns of 

(26) Qf = (Qi + Q2F)(I + F*F)-i/2 

(27) Ql = (Q2-QiF*)(I + FF*)-i/2 

are, respectively, orthornormal bases for invariant subspaces of . 
The representation of y^ with respect to QfjQI is, respectively, 

(28) T,l = (I + F*F)i/2p^^_r^^+^^2F](I + F*F)-i/2 

(29) S| = (I + FF*)-i/2[f22-Ffi2](I + FF*)V2 

where Si = diag{ai,a2, ■..,o'mm)- 

Corollary 1. Let p^, G (1/5, 1/4) be the only real root of the cubic polynomial p{p) = 1 — 8p + 
20p^ — 20p^ and .suppose 

(30) ^ < p,. 

Then Ran(Qf ) is the singular subspace associated with the s largest singular values of and 
Ran(Q2) the singular subspace associated with the rest of the singular values. 



Proof. It suffices to show that under (30) the smallest singular value of Si is larger than the largest 
singular value of S2. 

Since < 1/4, condition (30) implies that 

i|g||^„ „ < 1 

fmin — 2||<5||2 2 

which in turn implies (24). Note that under this condition ||F||2 < 1. By Proposition [2] we have 
(31) IISIII2 < ||I + FF*||2/^||£:22 -F£:i2||2 



< 1 + 7 " ' Il'?ll2 + 



(0"min - 2||<S||2)^y V 0"min-2||<S 



On the other hand, 

(32) min ||Sfe||2 = min ||(Si + + £:i2F)e|| ||I + F*F||2 

||e||2=l l|e||2=l 

> h7mm-h?2 " 1 + 



min 2||£^||2/ \ (fmin 2||£^||2)^ 



Let 

ll-^lb 



Imposing that the right hand side of (32) is greater than that of (31) leads to the inequahty 

p{p) = (1 - 2p)3 - 2p(l - 2p)2 - V 

= 1 - 8/9 + 20p2 - 20p^ > 
which holds for p < p.^, the only real root of It is readily verified that p* € (1/5, 1/4). 



□ 



In view of Corollary[T| it is natural to call Ran(Q2) the noise subspace and Ran(Q|) the signal 
(or object) subspace. 

Let be the orthogonal projection onto the noise subspace and define the MUSIC imaging 
function for the noisy data 

r{v) ^ 



We are ready to state the first main result of the paper. 



Theorem 1. Let = y + 8 where y and £ are given by (21)-(22). 
Suppose (f)j. G Ran($) if and only ifrGS. Then the condition 

(33) MSl < A 1 1 ' 



O-min 2 2 y^^r^ + 1 

implies that the s highest peaks of J^(r) coincide with the true locations of objects. Indeed, the 
object locations can be identified by the thresholding rule: 

(34) {r G /C : J%r) > IT'^} . 

Proof. Let {vj + 6vj : j = s + 1, n} be the columns of Ql. Clearly, 6vj is the (j — s)-th column 
of Ql — Q2- Now we have 

IIQI-Q2II2 < ||QiF*(I + FF*)-^/2||^^||Q^(j^pp*)-i/2_Q^||^ 
< ||QiF*||2 + ||(I + FF*)-i/2_j||2 

whose first term is bounded by ||F||2 and whose second term is bounded by 

||(I + FF*)-i/2-I||2 = ||(I + FF*)-i/2(j^(j^pp*)i/2)-ipp*||^ 

Hence 

(35) ||<5v,-||2 < ||Q|-Q2||2< ||F||2 + ^||F||2< ^^il^, j = s + l,...,n, 
where p = ||<?||2/o-min- 



By Corollary [T] { Vj + Svj : j = s + 1, ...,n} is the set of singular vectors associated with the n — s 
smallest singular values of Y. By definition, 



2 



k=s+l 
n 



(36) = E l^fc<^r|' + 2 E 3^[v^.0r</';5vfe] + ||(Q|-Q2; 

k=s+l k=s+l 



By assumption the first two terms on the right hand side of (36) vanish if and only if r G S. By 



(35) the third term is bounded by 



||(Q2-Q2) M2 < (i_2^)4 • 

For r 5, 



E |K+K)0r|' = E |v^'^r|'-||(Q|-Q2r</'ri 

fc=s+l 



(1 - 2p)4 



by (35). Hence J^{v) has the following behavior 



^2/1 ^n2\ -1 



(38) J'ir) < (n-'-^^^^z^) , reS^ 



Setting 



we obtain the inequality 



V(l-p)2 ^ ^ 5 (l_2/,)4 



(39) p'^-p+ ^^-^^ > 

4TS + 2V2 



whose solution is (33). Note that A < 1/5 < /j*, VF^ G [0, 1]. 



□ 



Condition (33) would not be very useful unless ||i?||2 can be bounded from above and fTniin, 
can be bounded from below by other known or accessible quantities. This is what the compressed 
sensing techniques enable us to do. 

3. Compressed sensing analysis 

We now give a quantitative evaluation of MUSIC based on compressed sensing theory. 
A fundamental notion in compressed sensing is the restrictive isometry property (RIP) due to 
Candes and Tao [3j. Precisely, let the sparsity s of a vector Z G be the number of nonzero 
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components of Z and define the restricted isometry constants G [0, f],(^+ G [0, oo) to be the 
smallest nonnegative numbers such that the inequahty 

(40) {i-6;)\\z\\l<\\^z\\l<{i + 6t)\\z\\l 

holds for all Z G C'^ of sparsity at most s. 

Roughly speaking this means that $ acts like a near isometry, up to a scaling, when restricted 
to s-sparse vectors. In particular, if d~_f_i < 1 then any s + 1 columns of $ are linearly independent 
which implies the characterization (13). 

More generally, let us extend the notion of the restricted isometry constants to ones 6^ associated 
with a particular set S, namely the smallest nonnegative numbers satisfying 

(41) {l-6s)\\Z\\l<\\^Z\\l<{l + 5+)\\Z\\l 

for all Z G supported on the set S. This will become important later when we analyze the case 
of an arbitrarily refined grid (Section [7]). Clearly, 

(42) (5^ = maxj|. 

\S\=s 

Then l\13n is equivalent to 6^, < 1 for all S' which is the union of S and another point r G S'^. 
First, let us estimate the magnitude of the error term £ in terms of E as follows. 



Lemma 1. Suppose (41) holds for ^. Then 



(43) ||£:||2 < ||E||2 + 2CmaxVl+4l|E||2 < ||E||i + 2Cmax Vl + 5.+ ||E||2 

where Cmax = ||Z||2 is the largest singular value ofZ. 

For the case of Z = X^'* with ^ satisfying (41), we have 

(44) ||£:||2 < ||E||2 + 2W(1 + 4)I|E||2 < ||E||2 + 2^^...^ + '^s+)l|E||2 
where ,^max = maxj 



Proof. First we have 



The RIP (41) then implies that 



\£h < l|E||i + 2||*ZE*||2. 



|*ZE*||^ < (l + 5+)||ZE*||^ 



and thus 



Sh < ||E||^ + 2a/1 + 5+||ZE*||2 



< ||E||^ + 2^1 + 5+||Z||2||E* 
In the case of scattering objects Z = X^* 



|ZE*||2 = ||E*X*||2 < ||E||2||*X*||2 < ||E||2W\/l+<^.+ 



s 



provided that ^ also satisfies the RIP (41). In this case. 



(45) ||^:||2 < ||E||2 + 2^n,ax(l + 5J)||E||2 



and hence (44). 



□ 
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Lemma 2. The minimum nonzero singular value fTmin of y obeys the lower hound 

(46) a^in > (1 - 5^ )Cmin > (1 " OCmin 

where 

II Z*e|| 2 

(47) Cmin = min ^ . 

eec= ||e||2 

For i/ie case of scattering objects Z = X^'* with ^ satisfying we have 

(48) a^in > (l-'55)'emin>(l-C)'din 

Proof. Using the max-min theorem [13] 

(49) cT,m = max min 

dimH=see'H ||e||2 

with % = Ran($) we obtain 

(50) C7mm(3^) > min ||*ZZ***e||2. 

eeRan(*) 

|e|l2 = l 

Let {uj : j = 1, s} be the eigen- vectors of associated with the nonzero eigenvalues {X\ > 
A2 > ... > Ag} and form the orthonormal basis of T-l. Write e = Ylj=i ^j^j where |ejp = 1. We 
have 

s 

and thus 

(51) ||**e||i = ^A,2|e,|2>A^ 



It follows from (50), (41 ) and (51 ) that 



(52) a^in >\/l-Ss\\ZZ*^*e\\2 >Vl-55CminA. 



by (47). On the other hand, A^ is exactly the smallest eigenvalue of G qsxs hence by 



(41) is bounded from below by 1 — 5^. Using this observation in (52) we obtain (46) 
In the case of scattering objects we can bound (^'^■^^ as 

Cmin — (l~^5)'^min 
Cmax — (-'- "I" "^5 )^max 



by using (41) with ^ and hence the result (48) 



□ 

Next we derive a lower bound for Ts in terms of RIC. 
Lemma 3. Fix S,\S\ = s. Then the lower bound is valid 

(53) > 1 — max — — ^ — — > 1 ^ • 
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Proof. Without loss of generality, suppose S = {l,2,...,s} and consider (pr = ^s+i- Let S' = 
5 U {s + 1}. Our subsequent analysis is independent of these choices modulo inconsequential 
notational change. 

Denote V(j)r = <p' and write the orthogonal decomposition 



(54) 

Hence we can express (j)' as 



b' = ^s+i - Cj^j = ^Z, Z = (-C1, -C2, -Cs, 1, 0, 0)^ G 



Using (41) for sparsity S' we obtain a lower bound for \\(j)'\\2 
(55) 



(l-S-,){l + Y,\cj\')<\\<P'\\l 
i=i 

On the other hand, we have by the Pythagorean theorem that 



(56) 



Applying (41) for sparsity 5 to the second term on the right hand side of (56) we obtain 

|2 



li-ll'A'lli = llE^A-ll2<(i + 4)Elc^i 



and hence 
(57) 



^|c,p>(l + 4)-i 
i=i 



||2 
r|l2 



Combining (57) and (55) we obtain 

/,'l|2 



4 



which can be solved to yield 
(58) 



A'l|2 



> 1 



2 + 6+ -5 



Minimizing (58) over r G 5"^ we obtain the first inequality in (53). 



The second inequality (53) follows from ( |42[ ) and the observation (by differentiation) that the 
quantity 

2 + 6+- 6^, 

is an increasing function of 6^, G [0,1] and 6^ G [0, oo) separately. 

□ 

Combining the preceding results we have the following stability criterion for exact recovery by 
MUSIC. 
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(59) 



If the noise-to-object ratio (NOR) satisfies 
e 



Theorem 2. Suppose 6^^^ < 1 (implying (13)) and ||E||2 = e. 

Cmax /Z \ 7+ 
-V 1 + oJ 



Cn 



C2 



+ {l-6s)A 



Cn 



where A is given by (33) then the object support S can be identified by the thresholding rule 



(60) 



r G /C : J^(r) > 2 1 



2 + 5t - 6] 



+1 



In the case of scattering objects Z = X^* with satisfying the RIP (40) the thresholding rule 



60) holds under the following bound on the noise-to-scatterer ratio (NSR) 

'^max 



(61) 



< J(l + 5+)2^ + (1 - 57)' A - (1 + 5J 

?min V ?min 



where ^max/Cmin is the dynamic range of scatterers. 



Proof By (43) and (46) 



(62) 



p' + 2pp^VT+St <A, 

Cmin 



Cn 



implies (33) in Theorem [T] The sufficiency of (59) now follows from solving the quadratic inequality 
(611) for p. 



The derivation of the thresholding rule (60) under the stronger condition (59) is exactly the same 



as that of (34). Alternatively, we can use (37) and (38) to verify validity of the thresholding rule 



(60) as follows. Let 



L = 1 



2 + 6t-S 



+1 



By using ( 38 ) and ( 39 ) it is straightforward to check that 2L ^ is greater than the right hand side 
of (38 1. On the other hand, (33) and Lemma [s] imply that 2L~' is smaller than the right hand side 
of fSTp. 

□ 



'he proof for the case of scattering objects is exactly the same as above. 



Remark 3. The right hand side of (61) decreases as the ratio 
(63) i^-^^s?^ 

decreases. In the underresolved case (Section^, 5^ is close to 1, making the ratio (63) a small 
number. For a noise-to-scatterer ratio smaller than (63) the s scatterers can be perfectly localized 
by the MUSIC algorithm with thresholding. This is the superresolution effect. 

A simple upper bound for the RIC can be given in terms of the notion of coherence parameter 
/i($) defined as 



max 



Zl\'^H\'El\'^lj\' 
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Namely, ^($) is the maximum of cosines of angles between any two columns. 
The proof of the following well known result is elementary and instructive. 

Proposition 3. For any r G N, we have 

5±</i(^)(r-l). 

Proof. Calculating the quantity ||^Z||2 



ZII2, we have 



\Z\ 



Using the quadratic inequality 2ab < + lP' we obtain 



Therefore (40) is satisfied with 5^ < /i($)(r — 1). 



□ 



Remark 4. For s = 2 it follows from Proposition^ that 

6f < ^(*). 

Since fi is almost surely less than unity for randomly selected sampling directions, the MUSIC 
algorithm will find the true location of object in the absence of noise, if there is only one object. 

The coherence bound for the most general setting of random sampling directions is this. 

Proposition 4. [9J Suppose any two points in K, are separated by at least ^ > 0. Let Sk,k = 1, n 
be independently drawn from the distribution f^ on the {d — \) -dimensional sphere independently 
and identically. Suppose 



(64) 



for any positive constants a,K. Then $ satisfies the coherence bound 

/i(*) < X + 



n 



with probability greater than (1 — a)^ where satisfies the bound 

(65) <ct{l+uji)-'^^f%,oo, d = 2 

(66) < ci(l + a;^)-i/ii,oo, d = 3 

where \\ ■ \\t,oo is the Holder norm of order t > 1/2 and the constant ct depends only on t. 



Remark 5. Replacing s^ and f^ in Proposition^^ by ^'jdfc and f^, respectively, we have the 
same conclusion about ^. 



The constraint (64) on the number of search points in the computation grid /C is relatively weak 
and allows an extremely refined grid. However, to have a small x^ ^ can not be small compared to 
the wavelength. 

Suppose u)i > C'^n for d = 2 or coi > C^/n for d = 3 where 

Cl 



(67) 



C > — = — max 
yj2K 



{llf lll.oo, 
15 



1 1,00 



Then, according to Propositions |4j [3] and Remark [5j with high probabihty, for continuously differ- 
entiable distributions /^,/' we have 

(68) <5^ < Jf+i < 2y/2Ks/./n. 

Theorem [2] then imphes the following. 



Corollary 2. Suppose uj£ > C^n for d = 2 or uj£ > C^/n for d = 3 with C given by (61). 
Suppose that ^/n/s > 4:\/2K (hence 6f,6f_^^ < 1/2 by (68)) and that the NSR obeys 



Smm \ V ^min <,mm I 

where A is given in (33). Then under the assumptions of Proposition^the MUSIC algorithm with 
the thresholding rule 



(70) <{ r G /C : J"(r) > ^'^^ 



25 } 

recovers exactly the locations of s scatterers with probability at least (1 — a)^. 
The value 128/25 is arrived from the fact that 

max T z — S 0/8. 

'5.+ ,<5J+i<l/2 2 + (5j-(5,+l 

4. Planar objects: optimal recovery 

Let us consider the favorable imaging geometry where all the scatterers lie on the transverse 
plane z = 0. Furthermore, we consider the idealized situation where the locations of the scatterers 
are a subset 5 of a finite square lattice /C of spacing i 

(71) JC = {rj-.j = l,...,N} = {{pi£,P2e,0):puP2 = 1,...,Vn} , j = (pi - 1)Vn + p2. 

Hence the total number of grid points is a perfect square. 
Suppose we choose the frequency such that 

(72) uji = V2tt. 

Let a^, = {ikiflkj^k = l,...,n be independently and uniformly distributed random variables in 
[—1, 1]^ and set 

(73) Sfe = -^(afc,V2-|afc|2). 

Let the incident directions d/, / = 1, m be selected the same way but independently from Sk,k = 
1, ...,n. It can be proved that with m > s the corresponding sensing matrix ^ has rank s with 
probability one . 

With (72)-(73) and j = {pi — l)\/iV + p2 the scattering amplitude ^ for linear extended objects 
yields the following extended sensing matrix 

(74) = e-^^^^-P G C"^^ 

The matrix (74) is often referred to as the random partial Fourier matrix in compressed sensing 
theory. 

The following is a standard result about random partial Fourier matrix [18j. 
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Proposition 5. [18] Suppose aj,j = 1, n are independently and uniformly distributed in [—1, 1]*^, d > 
1. // 

(75) > C6-^s In^ s In In - 

mn 7 

for 7 G (0, 1) and some absolute constant C , then with probability at least 1 — 7 the random partial 
Fourier matrix defined by (74) satisfies the RIC bound 

(76) 5f < 5^. 

Remark 6. The result holds true for sampling points aj which are i.i.d. uniform r.v.s in the 
discrete set 



instead of [—1, 1]*^ where N^/'^ is assumed to be an integer. 

Assume for simplicity the plane wave incidence as before. Choosing (5* = 1/2 for sparsity s + 1 
in Proposition [5] and using Theorem [2] we obtain the following result. 



Corollary 3. Suppose that (72) is true and that = d^, A; = 1, ...,n with 

(78) >4(:7(s + l)ln2(s + l)lniVln-. 

mn 7 



If the NSR obeys (69) the MUSIC algorithm with the thresholding rule (70) recovers exactly the 
locations of s scatterers with probability at least 1 — 7. 

4.1. Paraxial regime. Here we would like to extend the scattering problem to the paraxial regime 
for the preceding set-up where, instead of n sampling directions, n point sensors located on the 
transverse plane z = zq measure the scattered field. 

We shall make the paraxial approximation for the Green function between the object plane z = 
and the sensor plane z = zq: 

/'79^ Ql^ y\ — +y^)/{'izo)^-iu)xi/zo^-iu)yri/zo^iuj{£/-+ri^)/(2zo) 

47rzo 

with s = (,^, Tj, Zq)) r = (x, y, 0). Denote 

Let Sfc = (^fc, 7]}^, ZQ),k = 1, ...,n he the locations of the transceivers. We have ^ = ^ where the 
extended matrix $ is given by 

(80) l'fc,z = (G(si,r,),--- ,G(s„,r,))'^, j = l,...,iV 



where rj G K, defined in (71). After proper normalization the extended sensing matrix $ can be 
written as the product of three matrices 

(81) ^ = D1AD2 

where 

Di = diag(e^'^(«^'+''^')/(2"«)) G C"^'^, D2 = diag(e^"(^'+^')/(2^°)) G C^^^ 
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are unitary and 



■-•nxN 



Now suppose that 
and write {ijirjj) 
where p G Z^. Then with 

(82) 

A takes the form 
(83) 



l,...,n are independently and uniformly distributed in ^ A A^2 



a.j ■ A/2 with aj G 



-1,1] ,j = l,...,n. Write also 
Xzq 



2 ' 2 J 

ixi,yi,0) = ipe,o) 



— TTia-p] 



G C 



nxN 



which is exactly the random partial Fourier matrix given in (74). Here and below A denotes the 
wavelength. 

Since both Di and D2 are unitary and diagonal, they leave both the ^^-norm and the sparsity 
of a vector unchanged. Therefore Proposition [s] and Remarkjojare applicable to ^ given in (81). 
Analogous to Corollary |3] we have 

Corollary 4. Suppose that [8^ is true and that there are n transceivers satisfying (78). If the 



NSR obeys (69), then the MUSIC algorithm with the thresholding rule (70) recovers exactly the 
locations of s scatterers with probability at least 1 — 7. 

In the paraxial setting (82) in Corollary |4] replaces the condition (72) in Corollary [sj Condition 



(82) is exactly the classical Rayleigh criterion for resolution which is, in this case, the grid spacing 
i. 



5. Comparison with basis pursuit 
In the standard compressed sensing theory, one usually considers the following data model 

(84) = ^Z + E£C'', ZeC^, \\E\\2<e 

where the data and the object are vectors and employs the relaxed minimization principle called 
the Basis Pursuit Denoising (BPDN) [U 0] 

(85) min llZ'lli, s.t. - ^Z'lb < e 



for reconstruction. The noiseless version e = of (85) is called the Basis Pursuit (BP). Note that 
BPDN uses only one column of the MUSIC model (20). 

When Z is not s-sparse, consider the best s-sparse approximation Z^^") of Z. Clearly, Z^^^ = Z 
if Z is s-sparse. 

Denote the BPDN minimizer by Z. When does Z give a good approximation to the true Z? 
Again, the RIP (40) gives a useful characterization [2]. 

Theorem 3. Suppose the RIC of ^ £ C"^^ satisfies the inequality 
(86) + + 1) '^2-. < 1 
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Then the BPDN minimizer Z is unique and satisfies the error bound 
(87) \\Z-Z\\2 < Cxs-^I'^\\Z - Z^%x^C2e 

where 



Ci 



Co 



2 + (^/2 - 2)^2, + 



2s 



2s 



Remark 7. The real-valued version with 6s = 5f of Theorem^is proved in [2\. The proof for the 
complex-valued setting follows the same line of reasoning with minor modifications. For reader's 
convenience and for the purpose of showing where adjustments are needed, the full proof for the 
complex-valued setting is given in the Appendix\B{ 

Remark 8. Theorem^ does not guarantee exact recovery of support when E 0. 

An alternative approach to BPDN is greedy algorithms such as the orthogonal matching pursuit 
(OMP). The exact recovery of support by OMP is established for any s-sparse object vector Z such 
that the noise-to-object ratio satisfies 

e 1 ,^sA 

< o+M*)(; 



Zmin 2 2 



where Zja\n is the smallest absolute nonzero component of Z [8j . In order for the right hand side 



of (88) to be positive, it is necessary that 



11 

(89 s<- + ^. 

Remark 9. In the case of planar objects, under the assumptions of Proposition^ (with 5.^ = \/2—l), 
BP yields the exact solution Z = Z for n = 0{s) sampling directions (or sensors) and just one 
incident wave, modulo logarithmic factors. In comparison, the performance guarantee for MUSIC 
in Corollary^ assumes n = 0{s) sampling and incident directions. 

Using Proposition [3] and Theorem |3j we obtain 

Corollary 5. // 

1 V2-I 



s < 



2 2/i(*) 



cf. (89), then (81) holds true. 



Remark 10. Under the assumptions of Proposition with continuously differentiable f^ and 
BP recovers the s— sparse object exactly Z = Z in the noiseless case e = for n = O^s^) and 
sufficiently high frequency. 

This is similar to the performance guarantee for MUSIC in Corollary However, the perfor- 
mance guarantee for MUSIC assumes 0{s'^) incident waves while the performance guarantee for 
BP assumes only one incident wave. 
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6. Spectral estimation and source localization 



Let us turn to the original application where the MUSIC algorithm arises, namely the source 
localization and the frequency estimation for multiple random signals. The two applications share 
almost exactly the same mathematical formulation. 

Suppose the random signal x{t) consists of random linear combinations of s time-harmonic 
components from the set 

{e-''-^^':u, = j^,j = l,...,N}. 

Let us write 

N 

(90) = 

i=i 

and assume that there is a fixed set S (i.e. deterministic support) of s nonzero amplitudes and the 
elements in the complementary set are zero almost surely. 
Consider the noisy signal model 

(91) y{t) = x{t) + e{t) 

where e(t) is the Gaussian white-noise. The question is to find out which s components are non-zero 
by sampling y{t). 

Consider random sampling times t^, k = 1, n which are i.i.d. uniform r.v.s in the set {1, A^}. 
Write Y = {y{tk)) G C",^ = (e(tfc)) G C" and Z = (aj) G C^. Then by ^ we have 

(92) = ^Z + E 

$ . = _]_^-^2^tuj/N ^^nxN 

cf . ( 84 ) . From the one-dimensional setting of Proposition [5] and Remark ([g]) we know that if ( 75 ) 
is satisfied with 5^^ = \pl — 1 then the RIC of ^ obeys the bound (76) with probability at least 
1 — 7. Applying BPDN to (92) we obtain the error bound ( |87[ ) with 0(s) data, modulo logarithmic 
factors. 

How does MUSIC perform in this case? The standard MUSIC proceeds as follows. Let Ry = 
E[yy*] G C"^", Viz = E[ZZ*] G C^^^ and = E[EE*] be the covariance matrices of Y, Z and 
E, respectively. Note that R^ is sparse and has at most rank s. 

Suppose the noise and the signal are independent of each other. Then we have 

(93) Ry - Rb = ^Rz** 



which is of the form (18). 

Assume that Hz has rank s. This is true, for example when are zero- mean, independent 
random variables. In this case, Hz = diag(E|aip) has exactly s nonvanishing diagonal elements. 

Let X G C*^* be the 5x5 submatrix of R^. By assumption, X has rank s. Let $ = ^ = 
$5 G C"^'^ be the column submatrix restricted to the index 5, the (deterministic) support of Z. 
Let Y = Ry - Rs G We then can rewrite ^ in the form ^ suitable for MUSIC. 

Since X has full rank, the the ranges of $ and Y coincide. To guarantee the exact recovery 
by MUSIC, it suffices to show that the RIC of $ satisfies the bound < 1 which follows from 
Proposition ([s]) under the condition n = 0{s), modulo logarithmic factors, cf. (78). 

Let us formally state the performance guarantee for MUSIC as applied to the problem of spectral 
estimation. 
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Corollary 6. Suppose Rz has rank s and suppose that the noise is independent from the signal. 
Assume that the number n of time samples satisfies 

n 



(94) 



Inn 



> C{s + l)\n^ (s + l)lniVln7-\ 



cf. (78). Then for any noise level the singularities of the MUSIC imaging function (14) coincide 



with the frequencies in the random signals (90) with probability 1 — 7. 



Remark 11. The number of time samples assumed here is similar for both MUSIC and BPDN. 
However, many realizations of Y and Z are needed to calculate the covariance matrices accurately 



and form the equation (93) before the MUSIC reconstruction. Once (93) holds with sufficient 



accuracy, then the noise structure does not affect reconstruction as long as the noise is independent 
of the signal. 

In the absence of abundant realizations of signals, though, BPDN is the preferred method for 
spectral estimation. Indeed, BPDN can identify the frequencies approximately with just one re- 
alization of signals. The recovery error is at worst linearly proportional to the noise level as in 



(81). 



The source localization problem can be treated in the same vein as follows. 

Let us assume that s source points are distributed in the grid K. defined in (71) and each source 



point emits a signal described by the paraxial Green function ( 79 ) times the source amplitude which 
is recorded by the n sensors located at a^, z = 1, n in the plane z = zq. 

Let Z = (^(rj)) S C^. After proper normalization, the data vector Y can be written as (92) 
with the sensing matrix $ of the form (81). 



By the same analysis as before we arrive at the conclusion 
Corollary 7. Suppose has rank s and suppose that the noise is independent of the signals. 



Let the number n of time samples satisfy (94). For any noise level the singularities of the MUSIC 
imaging function (I4) coincide with the source locations with probability 1 — 7. 



7. Resolution and grid spacing 

Being essentially a gridless method, MUSIC 's flexibility with grid spacing is an advantage that 
the current BPDN-based imaging methods do not yet possess. 

Let £ a length scale to be determined below and let 5^ = {r G /C : dist(r,5) < £} be the i- 
neighborhood of the objects. For the problem of inverse scattering, dist(r,5) typically refers to the 
physical or Euclidean distance in the spatial domain. We would like to derive a thresholding rule 
which can eliminate all false alarms (i.e. artifacts) occurring outside Si, no matter how refined the 
grid spacing is relative to the frequency. 

Let ^ be the extension of $ over a fine grid of spacing i which may be much smaller than lo^^. 
When i = 0, the computation domain /C is a continuum. 



Generalizing the definition (19), we define 



(95) 



Ts{e) = mm \\cPr\\2^\\VcPr\\2 = /l-max||</.r||-20j:**t 



Clearly, Ts = Ts{0~^). In the noiseless case, the exact recovery of S by MUSIC is equivalent to 
Ts{0+) < 1. 

Extending the proof of Lemma [3] to {£) we have 
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Lemma 4. 

5^,(1 + 4; 



(96) Ts{£) > 1 - max 



S' = SeU{r} 2 + Jo - 6t, 

res? 



Extending the analysis leading to ([37|)-(p8|) we have 



V(i-p )^ 

(1 - 2p) 

The following result is analogous to Theorem [T] 



(98) J%r) < (riffl- ] , vest. 



Theorem 4. Suppose 5g, < 1,5' = 5 U {r}, Vr G SI 
If 

(99) iLM<A, = --- — 

(7rx\\r\ ^ ^ 



_ V2r5(£) + i 

then 

(100) 5 C e = {r : J"(r) > 2r^2(£)} 

w/iere 9 n = 0. 

Lemmas |4j [2] and Theorem [4] then implies the following result analogous to Theorem [2] 
Theorem 5. Suppose 5^, < 1,S' = S U {r}, Vr G Sf. If the NSR obeying the upper bound 



^2 f 

)Smax I /i r— /i i r+N^- 



max 



(101) _^ < (1 + 4)2^ + (1 - 5-)2A, - (1 + 5^ 

?min y Smin ?min 

TOt/j Ai as in Theorem then 



(102) 5 C G = < r : J^(r) > 2 1 - max 

w/iere 9 n 5f = 0. 




Let us now give an estimate of the length scale £ for (101) to be a useful upper bound for NSR. 
Let us focus on the general setting of Proposition [4j namely arbitrarily located scatterers and 
random sampling directions. 

We resort to the following result analogous to Proposition [3j The proof is exactly the same as 
before and is omitted here. 

Proposition 6. For any set B C IC, \B\ < r, we have 

4</x(*i5)(r-l). 



To proceed, let us tailor the estimate in Proposition |4] to the current setting as follows. 
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Proposition 7. [9j Suppose the physical distances between two points corresponding to any two 
members of B <Z K. are at least t. Let s^, A; = 1, ...,n be independently drawn from the distribution 
f^ on the {d — 1)- dimensional sphere independently and identically. Suppose 

I I - g 

for any positive constants a,K. Then satisfies the coherence bound 



n 



with probability greater than (1 — of' where satisfies the bound (65)-(66). 

Suppose w£ > C^n for d = 2 oi ujI > C^fn for d = 3 where C is given by (67) and assume that 
the s scatterers are separated by at least I from one another. Then, according to Proposition [7] and 
Proposition |6] with high probability, for any continuously differentiable sampling distribution 

b% < 6^, < 1^/2Ks|^/^ 

for all 5' = 5 U {r}, r G 5^ . Hence we have the following analogous result to Corollary^ 



Corollary 8. Suppose ujI > C^n for d = 2 or uj£ > C\pn for d = 3 with C given by l6''i_ 

Under the assumptions of Proposition^ (for B = S,S' ), y/njs > Ay/2K and the NSR bound 



§3' 

5 C e = <! r G /C : riv) > 
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with probability at least (1 — a)^ where B H 5^ = 0. 

8. Numerical tests 

In the simulations, zq = 10000, A = 0.1 and the search domain is [—250, 250]^ with grid spacing 
= 10 on the transverse plane z = 0. The scatterers are independently and uniformly distributed 
on the grid with amplitudes independently and uniformly distributed in the range [1,2]. The 
sensors are independently and uniformly distributed in the domain [—^4/2,^/2]^ with various A. 



The source locations are identical to the sensor locations. In the set-up, condition (82) is satisfied 
with A = 100. With these parameters, the paraxial regime is about to set in (cf. [llj). Note, 
however, that all the simulations are performed with the exact Green function. 

In our simulations we have used the Matlab codes YALLl (acronym for Your ALgorithms for 
LI, available at http://www.caam.rice.edu/ optimization/Ll/Y ALLl/| ). YALLl is a Ll- 
minimization solver based on the Alternating Direction Method [22]. 

Figure [2] compares the performances of MUSIC and BP in the well-resolved case ^ = 100 and 



the under-resolved case A = 10 where the aperture is only one tenth of that satisfying (82). For 
Figure [2] BP is carried out on the data matrix Y with the sensors coincident with the sources, 
i.e. $ = To put the problem in the proper set-up for BP, we vectorize Y by staking its n 
columns and denote the resulting vector by Y . We vectorize the diagonal matrix X by listing 
its iV diagonals as a C''^ vector X. The BP performance of this set-up has been analyzed in [llj. 
The numbers of recoverable scatterers shown in Figure [2] are computed for at least 90% recovery 
rate based on 100 independent realizations of transceivers and scatterers. In both cases, MUSIC 
recovers s = n — 1 scatterers with certainty. Clearly, for the well-resolved case, BP has a far 
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Figure 2. Comparison of MUSIC and BP performances, with both using the whole 
data matrix: the number s of recoverable scatterers versus the number of sensors 
n with A = 100 (left), the well-resolved case, and A = W (right), the under- 
resolved case. In the well-resolved case, BP delivers a much better (quadratic-in- 
77.) performance than MUSIC; in the under- resolved case, MUSIC outperforms BP 
whose performance tends to be unstable in this regime. The numbers of recoverable 
scatterers by BP are calculated based on successful recovery of at least 90 out of 
100 independent realizations of transceivers and scatterers while the success rate of 
MUSIC is 100%. 
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BP 
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Figure 3. Comparison of MUSIC and BP performances with BP employing only 
single column of the data matrix: the number s of recoverable scatterers versus the 
number of sensors n with A = 100 for n G [10,30] (left) and n £ [150,200] (right). 
Both BP curves show a roughly linear behavior with slope less than that of the 
MUSIC curves. 



superior performance to MUSIC. Indeed, it can 
with high probability in the well-resolved case 
near-parabolic curve in Figure Islleft). For the 



be shown that BP can recover s = 0{v?) scatterers 
[11]. The quadratic behavior is illustrated by the 
under-resolved case, however, MUSIC outperforms 
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Figure 4. Success probability of the MUSIC reconstruction versus aperture for 
n = 10,s = 9 (left), n = 100, s = 9 (middle) and n = 100, s = 99 (right). Note 
the different aperture ranges for the three plots. The success rate is calculated from 
1000 trials. Increasing the number of transceivers for the same number of scatterers 
reduces the aperture required for the same success rate. The reduction of aperture 
is about three folds (left to middle). On the other hand, higher number of scatterers 
with the same number of transceivers also demands larger aperture for the same 
success rate. The increase in aperture is about 7 times (middle to right). 
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Figure 5. Success probability of MUSIC versus the number of transceivers with 
A = 0.5, s = 9 (left), A = 0.2, s = 9 (middle) and ^ = 15,s = 99 (right). The 
probabilities are calculated from 1000 independent trials. 



BP by a significant margin. Figure [sj^right). As pointed out in Remark [sj the MUSIC algorithm 
has the superresolution capability for a sufficiently small noise-to-scatterer ratio. 

If only one column of Y is used in BP as discussed in Section [5| then MUSIC outperforms BP 
by a wide margin even in the well-resolved case. Figure [3] 

We further investigate the performance of the MUSIC algorithm for the extremely under-resolved 
case when BP essentially has extremely low probabilities of exact recovery (even for s = 1). Figure 
|4] shows the success probabilities of MUSIC as a function of aperture for various n and s while 
Figure [5] shows the success probabilities of MUSIC as a function of n for various A and s. The 
success rates are calculated from 1000 independent realizations of transceivers and scatterers. 

Three observations about Figure [4] are in order: (i) The optimal performance of s = n — 1 does 
not hold with certainty for s relatively large with respect to aperture (cf. left and right panels); (ii) 
Increasing the number of randomly selected transceivers reduces the aperture required for the same 
probability of recovering the same number of scatterers (left to middle panels); (iii) Increasing the 
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A=10 




Figure 6. Success probability of MUSIC reconstruction of s = 10 scatterers with 
n = 100 transceivers versus the noise level a in the well-resolved case A = 100 
(left) and the under-resolved case A = 10 (right). The success rate is calculated 
from 1000 trials. Note the different scales of a in the two plots. Noise sensitivity 
increases dramatically in the under-resolved case. 



number of randomly selected scatterers increases the aperture required for the same probability of 
recovery with the same number of transceivers (middle to right panels). 

Likewise, the success rates increase with the number of transceivers for any aperture and sparsity 
(Figure[5]). The most interesting plot in Figure[5]is the middle panel which shows for A = 0.2, s = 9 
the success rate curve becomes a plateau after reaching 80%. This is not inconsistent with the 
prediction of Proposition [l] since Proposition [l] assumes a fixed configuration of scatterers while 
Figure [5] is for random, independent realizations of scatterers. In other words the threshold uq in 
Proposition [T] may not be uniformly valid for all configurations of s scatterers in the under-resolved 
case. On the other hand, when the aperture increases by two and half times to ^ = 0.5 and 
the number of transceivers increases to 15, the performance becomes uniform with respect to the 
scatterer configuration (left panel). 

Figure[6]shows the noise sensitivity of MUSIC reconstruction of 10 scatterers with 100 transceivers. 



Here n and s are chosen so that (75) is roughly satisfied. We add the i.i.d. noises 



(103) a{ei + ie2)ymax 

to the entries of the unperturbed data matrix where ei and 62 are independent, uniform r.v.s in 
[— 1, 1] and ymax is the maximum absolute value of the data entries. Hence the signal-to- noise 
ratio (SNR) is about 2~^a~'^. In the well-resolved case {A = 100) the MUSIC reconstruction can 
withstand a significant amount of noise in the data matrix. Indeed, at SNR 0.5 the success rate 
is almost 100%, consistent with the prediction of Theorem [2| and even at SNR 0.22 (cr = 1.5) the 
success rate can be indefinitely improved by increasing the number of transceivers (Figure [7| left 
panel) . 

In the under-resolved case, however, the noise sensitivity increases significantly. Figure [6] (right 
panel) reminds us how fragile the superior performance of MUSIC in the under-resolved case is, 
cf . Figure [s] (right panel) . Figure [T] (right panel) further indicates that in the under- resolved case 
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A=100,s=10,a=1.5 
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Figure 7. Success probability of MUSIC reconstruction of s = 10 scatterers as a 
function of n with a = 150% in the well-resolved case A = 100 (left) and a = 5% 
in the under-resolved case ^4 = 10 (right). The success rate reaches the plateau of 
85% near n = 1000 in the under-resolved case. The success rate is calculated from 
1000 trials. 



the success rate may not be indefinitely improved by increasing the number of transceivers in the 
presence of noise. 



9. Conclusion 

We have developed a framework for discrete, quantitative analysis of the MUSIC algorithm in 



the well- resolved case. Our approach is based on the RIP (40) and its variant (41 ) which takes into 
account of the object configuration as well as sparsity. 

Our first main result is a support recovery condition (Theorem [2]) that for the NOR obeying (59) 
MUSIC can exactly localize the objects with noisy data. Our result indicates the superresolution 
capability of the MUSIC algorithm when the noise level is sufficiently low (Remark ^ . 

We have provided a coherence approach to estimating RIC (Propositions [4] and |3|) for general 
object configuration in three dimensions with the grid spacing £ ~ As and the sensor number n ~ s^. 
When the scatterers are distributed in a transverse plane, then £ ~ A, n ~ s (modulo logarithmic 
factors), suffices. We have extended these results to the gridless setting for which i is interpreted 
as the minimum distance between objects and only approximate localization up to the error i is 
sought (Theorem [5] and Corollary [s]) . 

Our comparative analysis shows that when the whole data matrix is employed in both BP and 
MUSIC, BP outperforms MUSIC in the well- resolved case in the sense that the number of objects 
recoverable by BP grows quadratically with the number of transceivers while that by MUSIC grows 
linearly. The MUSIC reconstruction can tolerate a significant amount of noise in the data matrix 
(Figures |6| [7] left panels). On the other hand, our numerical results show that in the under-resolved 
case MUSIC outperforms BP by a wide margin (Figure [2| right panel, Figures [4] and [s]) . However, 
music's superresolution effect is still unstable with respect to noise in the data matrix (Figures 
[6|[7| right panels). 
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Finally, even in the well-resolved case where the employment of just one column of the data 
matrix by BP guarantees a probabilistic recovery of objects numbered in linear proportion to the 
number of sensors, analogous to the performance guarantee of MUSIC, the latter outperforms the 
former in numerical simulations by a wide margin (Figure [s]). 
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Figure 8. Scattering by extended objects 

Appendix A. Sparse extended objects 

In this appendix we extend the MUSIC algorithm to image sparse extended scatterers by inter- 
polating from grid points. 

Suppose that the object function ^(r) has a compact support. Consider the discrete approxima- 
tion by interpolating from the grid points 

eKr)=^' j;5(r/^-qm), X C Z'^ 

where g is some spline function and £ is the grid spacing. Since has a compact support, X is a 
finite set. For simplicity assume d = 2 and let X be the finite lattice 

2^ = {q= {(11,(12) ■■ qi,q2 = 1,:.,Vn} 

of total cardinality and IC = £1. In the case of a characteristic function g, is a piece- wise 
constant object function. We will neglect the discretization error and assume ^^(r) = ^(r) in the 
subsequent analysis. 

The data matrix Y G ([^i^x'm jg giyg^ by 



(104) Yk,i ~ £^Ta£ci) f 5(r7^-q)e^"('^'-'^)"-'dr', 



= £''{27rf^'g{£u{di - s,)) ^ C{i^)e''''^^'-'^^-^. 

As before we maintain the option of normalizing Y. Suppose {^£ = ?(^qj) '■ j = 1, •••) is the set 
of nonvanishing ^(^q) and let X = diag(^j) G C'"'. Dividing jioit by £'^{27r)'^/'^g{£uj{di - Sfe))/2 
we can write the data matrix in the form ( |10[ ) with the sensing matrices 



n 
1 



Jn 

where j = (gi - 1)\/A + q2 
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In other words, the scattering analysis for both point and extended scatterers leads to the same 
type of Fourier- like matrices. 

Appendix B. Proof of Theorem [3] 
The following lemma differs from the original version in p]. 
Lemma 5. We have 

1 



for all Z,Z' supported on disjoint subsets T,T' C {1, ...,m} with \S\ < s, \S'\ < s' . 

Proof. Without loss of generality, suppose Z,Z' are unit vectors. Since Z _L Z' , \\Z ± Z'\\\ = 2. 
Hence we have from the RIP (40) 

(105) 



2(1 - 6-.,) < mz ± Z')\\i < 2(1 + 5+ 



By the parallelogram identity and ( 105 ) 
which proves the lemma. 



^z + ^z'\\l-\\^z -^z'wl 



□ 



By the triangle inequality and the fact that Z is in the feasible set we have 
(106) ||*(Z - Z)\\2 < W^Z - Y\\2 + \\Y - *Z||2 < 2e. 

Set Z = Z + A and decompose A into a sum of vectors A^^, A^^, A^j, each of sparsity at most 
s. Here Sq corresponds to the locations of the s largest coefficients of Z; Si to the locations of the 
s largest coefficients of A^g; 52 to the locations of the next s largest coefficients of A^g, and so on. 
Step (i). For j > 2, 



lAcJIo <si/2||A.,|U < s-i/2||A, 



'j-i I 



and hence 
(107) 

This yields 
(108) 

Also we have 

ll^lli > ll^lli 
which implies 
(109) 



Y,\\As^h<s-y'Y.\\^sMi<s-'/'\\As^,h. 
i>2 i>i 



H5oUS'i)'=l|2 - II '^^Sjh < ll^5jl|2 < S ^/^IIAsglli. 

i>2 j>2 



\Zso + ^5olll + ll-^Sg + ^SljWl > l|-^5olll - II^Solll 
IIA^glli < 211^5^111 + IIAs'olli. 



I^5g||l 



I Ac 



Note that yZ^gUi = ||Z — Z^^'^Hi by definition. Applying (108), (109) and the Cauchy-Schwarz 
inequality gives 

(110) ||A(5„u5i)=ll2 < ||A5ol|2 + 2eo 

where cq = - ZW||i. 
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Step (ii). Observe 



(111) ||*A 



i>2 



= 5r(*Asou5i,*a) - J] [3?(*A5o,*A5,) +K(*A5i,*As, 

This calculation differ slightly from the corresponding calculation in [2]. 
From (fToel) and the RIP dtol) it follows that 



*AsoU5i/>|| < II^A^ouSilbll^Alla < 2eY'l + <5J,||A5ouSlll2• 
Moreover, it follows from Lemma [B] that 

k(*A5o,^A5,)| < l{6l + 6^^)\\AsM^s,\\2 

K(*A5„*A5,)| < ^(<^2';+'^2'.)l|A5ol|2||A5,||2 

for j > 2. Since Sq and are disjoint 



|Asol|2+ IIA5JI2 < V2J||AsJ|2+ IIA5JII = \/2||A5„u5il|2. 



Also 



(l-'^2"JI|A5oUSill2 < ll*A5oU5il 



< II Asou5, II2 (^e^Jl + 5+ + {51 + ) Yl 11^^ 



Therefore from (107) we obtain 



2Wl + (5+ 

l^SoUSih < + ps ^'^\\As^\\i, a = — 



l-(5, 



2s 



72 (^^2^^ + ^2-.) 

l-<5;" 



and moreover by ( 109 ) 
Namely, 



IA50US1II2 < ae + p||Asg||2 + 2/360. 
||Asou5ill2 < (1 - p)~Hae + 2y9eo) 



if (jSej) holds. 
Finally, 

IIAII2 < ||A5„u5j|2 + ||A(s„u5i)-ll2 < 2||AsoU5il|2 + 2eo < 2(1 - p)-\a£ + (1 + p)eo) 
which is what we set out to show. 
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